/*
  File:    WENO.h
  Purpose: weighted ENO (essentiall non-oscillatory) schemes for discretizing 
           |Speed * grad(u)|   

  Contact: Paul Macklin
           pmacklin@math.uci.edu
		  http://math.uci.edu/~pmacklin
*/

#ifndef WENO_H_
#define WENO_H_

// #include "../ImprovedMath.h"
#include "../Matrix.h"

extern double _EPS;
extern double dx;
extern double dy;

double WENOweight( double dA, double dB, double dC, double dD);
Matrix WENO3(Matrix &mSpeedIn, Matrix &mIn);
Matrix WENO5(Matrix &mSpeedIn, Matrix &mIn);
Matrix WENO5banded(Matrix &mSpeedIn, Matrix &mIn, double BandSize );
Matrix WENO5VectorSpeed(Matrix &mSpeedInX, Matrix &mSpeedInY, Matrix &mIn);

#ifndef DIFF_QUOTIENTS
#define DIFF_QUOTIENTS

inline double ForwardDiffX(Matrix &mIn, double dx, int nI, int nJ)
{ return ( *mIn(nI+1,nJ) - *mIn(nI,nJ) ) / dx; }

inline double BackwardDiffX(Matrix &mIn, double dx, int nI, int nJ)
{ return (*mIn(nI,nJ) - *mIn(nI-1,nJ) ) / dx; }

inline double ForwardDiffY(Matrix &mIn, double dy, int nI, int nJ)
{ return ( *mIn(nI,nJ+1) - *mIn(nI,nJ) ) / dy; }

inline double BackwardDiffY(Matrix &mIn, double dy, int nI, int nJ)
{ return ( *mIn(nI,nJ) - *mIn(nI,nJ-1) ) / dy; }

inline double CentDiffX(Matrix &mIn, double dx, int nI, int nJ)
{ return ( *mIn(nI+1,nJ) - *mIn(nI-1,nJ) ) / (2.0*dx); }

inline double CentDiffY(Matrix &mIn, double dy, int nI, int nJ)
{ return ( *mIn(nI,nJ+1) - *mIn(nI,nJ-1) ) / (2.0*dy); }

inline double CentDiffXX(Matrix &mIn, double dx, int nI, int nJ)
{ return ( *mIn(nI+1,nJ) -2.0*(*mIn(nI,nJ)) + *mIn(nI-1,nJ) ) / ( square(dx) ); }

inline double CentDiffYY(Matrix &mIn, double dy, int nI, int nJ)
{ return ( *mIn(nI,nJ+1) -2.0*(*mIn(nI,nJ)) + *mIn(nI,nJ-1) ) / ( square(dy) ); }
#endif

#ifndef DIFFERENTIALS
#define DIFFERENTIALS
inline double ForwardDiffX(Matrix &mIn, int nI, int nJ)
{ return ( *mIn(nI+1,nJ) - *mIn(nI,nJ) ); }

inline double BackwardDiffX(Matrix &mIn, int nI, int nJ)
{ return (*mIn(nI,nJ) - *mIn(nI-1,nJ) ); }

inline double ForwardDiffY(Matrix &mIn, int nI, int nJ)
{ return ( *mIn(nI,nJ+1) - *mIn(nI,nJ) ); }

inline double BackwardDiffY(Matrix &mIn, int nI, int nJ)
{ return ( *mIn(nI,nJ) - *mIn(nI,nJ-1) ); }

inline double CentDiffX(Matrix &mIn, int nI, int nJ)
{ return ( *mIn(nI+1,nJ) - *mIn(nI-1,nJ) ) ; }

inline double CentDiffY(Matrix &mIn, int nI, int nJ)
{ return ( *mIn(nI,nJ+1) - *mIn(nI,nJ-1) ); }

inline double CentDiffXX(Matrix &mIn, int nI, int nJ)
{ return ( *mIn(nI+1,nJ) -2.0*(*mIn(nI,nJ)) + *mIn(nI-1,nJ) ) ; }

inline double CentDiffYY(Matrix &mIn, int nI, int nJ)
{ return ( *mIn(nI,nJ+1) -2.0*(*mIn(nI,nJ)) + *mIn(nI,nJ-1) ); }
#endif

#ifndef WENO_FUNCTIONS
#define WENO_FUNCTIONS 

inline double PlusFunction( double dIn)
{ if( dIn > _EPS ){ return dIn; }
  else{ return 0.0; }
}
inline double MinusFunction( double dIn)
{ if( dIn < -_EPS ){ return dIn; }
  else{ return 0.0; }
}

inline double IS0( double dA, double dB, double dC, double dD)
{ return 13.0 * square(dA-dB) + 3.0*square(dA-3.0*dB); }

inline double IS1( double dA, double dB, double dC, double dD)
{ return 13.0 * square(dB-dC) + 3.0*square(dB+dC); }

inline double IS2( double dA, double dB, double dC, double dD)
{ return 13.0 * square(dC-dD) + 3.0*square(3.0*dC-dD); }
#endif

#endif
